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Abstract. We study the origin of the bi-stability jump 
in the terminal velocity of the winds of supergiants near 
spectral type Bl. Observations show that here the ratio 
Voo/vosc drops steeply from about 2.6 at types earlier than 
Bl to a value of Voo/vcsc=l-3 at types later than B2. To 
this purpose, we have calculated wind models and mass- 
loss rates for early-type supergiants in a Tcff grid covering 
the range between T^s = 12 500 and 40 000 K. These 
models show the existence of a jump in mass loss around 
Tcff = 25 000 K for normal supergiants, with M increasing 
by about a factor five from Toff ~ 27 500 to 22 500 K 
for constant luminosity. The wind efficiency number rj ~ 
Mvaol {Lf/c) also increases drastically by a factor of 2 - 3 
near that temperature. 

We argue that the jump in mass loss is accompanied 
by a decrease of the ratio Voo/vosc, which is the observed 
bi-stability jump in terminal velocity. Using self-consistent 
models for two values of T^s, we have derived Vqo/vcsc = 
2.4 for Teff = 30 000 K and v^/vesc = 1-2 for T^s = 17 
500 K. This is within 10 percent of the observed values 
around the jump. 

Up to now, a theoretical explanation of the observed 
bi-stability jump was not yet provided by radiation driven 
wind theory. To understand the origin of the bi-stability 
jump, we have investigated the line acceleration for mod- 
els around the jump in detail. These models demonstrate 
that M increases around the bi-stability jump due to an 
increase in the line acceleration of Fe III below the sonic 
point. This shows that the mass-loss rate of B-type super- 
giants is very sensitive to the abundance and the ionization 
balance of iron. 

Furthermore, we show that the elements C, N and O 
are important line drivers in the supersonic part of the 
wind. The subsonic part of the wind is dominated by the 
line acceleration due to Fe. Therefore, CNO-processing is 
expected not to have a large impact on M but it might 
have impact on the terminal velocities. 
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Finally, we discuss the possible role of the bi-stability 
jump on the mass loss during typical variations of Lumi- 
nous Blue Variable stars. 

Key words: Radiative transfer - Stars: early-type - 
Stars: mass-loss ~ Stars: supergiants - Stars: winds 



1. Introduction 

In this paper wc investigate the origin and the conse- 
quences of the bi-stability jump of the stellar winds of 
early- type stars near spectral type Bl. This bi-stability 
jump is observed as a steep decrease in the terminal ve- 
locity of the winds from Voo — 2.6wcsc for supergiants of 
types earlier than Bl to Voo — l-3wcsc for supergiants of 
types later than Bl (Lamers et al. 1995). We will show 
that this jump in the wind velocity is accompanied by a 
jump in the mass-loss rate with M increasing by about a 
factor of five for supergiants with Toff between 27 500 and 
22 500 K. 

The theory of radiation driven winds predicts that the 
mass-loss rates and terminal velocities of the winds of 
early-type stars depend smoothly on the stellar param- 
eters, with Voo — 3tiosc and M cx L^-^ (Castor et al. 1975, 
Abbott 1982, Pauldrach et al. 1986, Kudritzki et al. 1989). 
This theory has not yet been applied to predict the ob- 
served jump in the ratio Voo/vosc for supergiants near spec- 
tral type Bl. The change from a fast to a slow wind is 
called the bi-stability jump. If the wind momentum Mv^o 
were about constant across the bi-stability jump, it would 
imply that the mass-loss rate would increase steeply by 
about a factor of two from stars with spectral types ear- 
lier than Bl to later than Bl. Unfortunately, there are 
no reliable mass-loss determinations from observations for 
stars later than spectral type Bl. 

So far, a physical explanation of the nature of this 
bi-stability jump has been lacking. In this paper, we at- 
tempt to provide such an explanation and we investigate 
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the change in mass-loss rate that is accompanied by the 
change in Voo- 

The concept of a bi-stabihty jump was first described 
by Pauldrach & Puis (1990) in connection to their model 
calculations of the wind of the Luminous Blue Variable 
(LBV) star P Cygni {Tes= 19.3 kK). Their models showed 
that small perturbations in the basic parameters of this 
star can either result in a wind with a relatively high mass 
loss, but low terminal velocity, or in a wind with relatively 
low M, but high v^c ■ Their suggestion was that the mecha- 
nism is related to the behaviour of the Lyman continuum. 
If the Lyman continuum exceeds a certain optical depth, 
then as a consequence, the ionization of the metals shifts 
to a lower stage. This causes a larger line acceleration 
and finally yields a jump in M. 

The models of Pauldrach & Puis (1990) for P Cygni 
show that the wind momentum loss per second, Mvoo, is 
about constant on both sides of the jump (see Lamers & 
Pauldrach 1991). So Lamers et al. (1995) put forward the 
idea that the mass-loss rate for normal stars could increase 
by about a factor of two, if v^o decreases by a factor of 
two, so that Mvao is constant on both sides of the jump. 

Whether this is indeed the case, is still unknown. To in- 
vestigate the behaviour of the mass loss at the bi-stability 
jump, we will derive mass-loss rates for a grid of wind mod- 
els over a range in Teff- The main goal of the paper is to 
understand the processes that cause the bi-stability jump. 
Although our results are based on complex numerical sim- 
ulations, we have attempted to provide a simple picture of 
the relevant physics. We focus on the observed bi-stability 
jump for normal supergiants. Nevertheless, these results 
may also provide valuable insight into the possible bi- 
stable winds of LBVs. 

It is worth mentioning that Lamers & Pauldrach 
(1991) and Lamers et al. (1999) suggested that the bi- 
stability mechanism may be responsible for the outflow- 
ing disks around rapidly-rotating B[e] stars. Therefore our 
results may also provide information about the formation 
of rotation induced bi-stable disks. 

The paper is organized in the following way. In Sect. 
^ we describe the basic stellar wind theory. In particular 
we concentrate on the question: "what determines M and 
Woo?". We show that M is determined by the radiative ac- 
celeration in the subsonic region. In Sect. ^ we explain the 
method that we use to calculate the radiative acceleration 
with a Monte Carlo technique and the mass-loss rates of 
a grid of stellar parameters. Sect, ^describes the proper- 
ties of the models for which we predict M. In Sect. || our 
predicted bi-stability jump in mass loss will be presented. 
Then, in Sect. ^ we discuss the origin of this jump and 
show that it is due to a shift in the ionization balance of 
Fe IV to Fe iii. Then, we discuss the possible role of the 
bi-stability jump in M on the variability of LBV stars in 
Sect. ^. Finally, in Sect. |^, the study will be summarized 
and discussed. 



2. What determines M and v^l 

2.1. The theory of M determination 

Mass loss from early-type stars is due to radiation pressure 
in lines and in the continuum (mainly by electron scatter- 
ing). Since the radiative acceleration by line processes is 
the dominant contributor, the winds are "line-driven" , i.e. 
the momentum of the radiation is transferred to the ions 
by line scattering or line absorption. Line-scattering and 
line absorption occur at all distances in the wind, from 
the photosphere up to distances of tens of stellar radii. So 
the radiative acceleration of the wind covers a large range 
in distance. 

The equation of motion of a stationary stellar wind is 
dv GM^. 1 dp , . ,^ , 

^^:r = — 2 :r+frad(?-) (i) 

dr p dr 

where (jrad is the radiative acceleration. Together with the 
mass continuity equation 

M = 47rr2p(r)u(r) (2) 

and the expression for the gas pressure p = TZpT / p,, where 
TZ is the gas constant and fi is the mean mass per free 
particle in units of mn, we find the equation of motion 

where a is the isothermal speed of sound. For simplicity 
we have assumed that the atmosphere is isothermal. In 
this expression the effective mass Meff = M,(l — Fe) is 
corrected for the radiation pressure by electron scattering. 
(7l is the line acceleration. The equation has a singularity 
at the point where v{r) = a, this critical point is the sonic 
point. If the line acceleration gh{r) is known as a function 
of r, the equation can be solved numerically. A smoothly 
accelerating wind solution requires that the numerator of 
Eq. ^ reaches zero exactly at the sonic point where the 
denominator vanishes. 

It should be stated that this critical point (sonic point) 
at Tc ~ 1.025i?* and Vc — 20 km s^^ is not the same as 
the CAK critical point. The CAK critical point is located 
much further out in the wind at rc ~ 1.5i?* and about 
Vc — O.bvoo- If the line acceleration gi^ in Eq. |3| were to 
be rewritten as a function of velocity gradient instead of 
radius, then one would find the CAK critical point. Paul- 
drach et al. (1986) showed that if the finite disk correction 
to the CAK theory is applied, then the Modified CAK crit- 
ical point moves inward and is located at rc — 1.04i?* and 
at Vc — 100 km s^^. This is much closer to the sonic point! 
Although the (Modified) CAK critical solution may well 
provide the correct mass-loss rate and terminal velocity, 
there is concern about its physical reality (see e.g. Lucy 
1998 and Lamers & Cassinelli 1999 for a thorough dis- 
cussion). Lucy (1998) has given arguments favouring the 



Jorick S. Vink et al.: The nature of the Bi-stability jump 



3 



sonic point as the physical more meaningful critical point. 
We will use the sonic point as the physically relevant crit- 
ical point. This is the point where the mass- loss rate is 
fixed. Throughout the paper we will therefore refer to the 
subsonic part of the wind for the region close to the pho- 
tosphere where the mass loss is determined, and to the su- 
personic part for the region beyond the sonic point where 
the mass-loss rate is already fixed, but the velocity has 
still to be determined. 

The critical solution can be found by numerically in- 
tegrating Eq. ^, starting from some lower boundary tq in 
the photosphere, with pre-specified values of To and po and 
with a trial value of vq. The value of vq that produces a ve- 
locity law that passes smoothly through the critical point 
is the correct value. Alternatively, for a non-isothermal 
wind with a pre-specified T(r)-relation, one can integrate 
inwards from the critical point with an assumed location 
Tc, and then adjust this value until the inward solution 
gives a density structure that reaches r — 2/3 at the lo- 
cation where T(r) = Tofi (e.g. see Pauldrach et al. 1986). 
The critical solution specifies the values of tq ~ i?*, pa 
(given by r(ro) — 2/3) and vq at the lower boundary. 
This fixes the value of M via the mass continuity equa- 
tion (Eq. H) . Note that M is determined by the conditions 
in the subsonic region! 

We will show below that an increase in gh(r) in the 
subsonic region results in an increase in M. This can be 
understood because in the subsonic region, where the de- 
nominator of Eq. 1^ is negative, an increase in gi^ gives a 
smaller velocity gradient. Integrating from the sonic point 
inwards to the lower boundary with a smaller velocity gra- 
dient, implies that the velocity the lower boundary should 
be higher and hence the mass-loss rate, M — AirrQPQVQ, 
must be higher. On the other hand, an increase in gi^ in 
the supersonic region, yields a larger velocity gradient and 
this would directly increase the terminal velocity Vao- 

Another way to understand how an increase in M is 
caused by an increase in gi^ below the sonic point, is based 
on the realization that the density structure of the sub- 
sonic region is approximately that of a static atmosphere. 
This can be seen in Eq. 0. Since the term v dv/dr is much 
smaller than the acceleration of gravity, it can approxi- 
mately be set to zero in the subsonic region. (This is not 
correct close to the sonic point.) In an isothermal static at- 
mosphere the density structure follows the pressure scale- 
height. Adding an extra outward force in the subsonic re- 
gion results in an increase of the pressure-scaleheight and 
hence in a slower outward decrease in density. This means 
that just below the sonic point, where w ~ a, the density p 
will be higher than without the extra force. Applying the 
mass continuity equation (Eq. ^ at the sonic point then 
shows that the mass-loss rate will be higher than with- 
out the extra force in the subsonic region. (See Lamers & 
Cassinelli 1999 for a thorough discussion). 




Fig. 1. Extra "bumps" on the radiative acceleration gi,{r) 
below the sonic point. The solid fine is ghir) of the model 
without a "bump". The dotted lines show gh{r) with the 
adopted bumps with peakheights of 150, 300 and 500 cm 
s~^. The cross indicates the sonic point at 1.0135 i?*. 

2.2. A simple numerical experiment: the sensitivity of M 
on the subsonic gL 

A simple numerical experiment serves to demonstrate the 
dependence of M on the radiative acceleration in the sub- 
sonic region. We start with an isothermal model of the 
wind from a star of M^s = 20Mq, = 16.92i?0, Tes 
= 25 000 K, Twind = 0.8Toff= 20 000 K. We then specify 
the line acceleration gh{r) in such a way that it produces 
a stellar wind with a mass-loss rate of 1.86 lO~^Af0yr~^ 
and with a /3-type wind velocity law 



v(r) 



(1 - 

r 



(4) 



where (3=1 and Voo = 1500 km s~^. (This ghir) is found 
by solving Eq. ^ with this fixed velocity law) . This model 
is very similar to one of the models near the bi-stability 
jump that we will calculate in detail in Sect. H. As a lower 
boundary we choose the point where p = 10 ° g cm~'^ at 
To = i?* . Figure ^ shows the resulting variation of (r) . 
Adopting this variation of gh{r) and solving the momen- 
tum equation with the condition that the solution goes 
smoothly through the sonic point, we retrieve the input 
mass-loss rate and input velocity law, as one would ex- 
pect. The sonic point is located at = 1.0135i?* where 
V = 16.6 km s~^, and where gi,{r) = 1.63 10^ cm s~^. 

Let us study what happens to M and v^o if we change 
the line acceleration in the subsonic region. To this pur- 
pose we add a Gaussian "bump" to gh{r). This bump is 
characterized by 



bump/ \ 



peak 

5l exp 




(5) 



where z = l/{(r/i?,) — 1}, Zp = 150 describes the location 
of the peak at r/R^ = 1.0067 and Az = 30 gives the width 
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Fig. 2. The effect of increasing tlie line acceleration in the 
subsonic region on M (upper panel) and a simple deriva- 
tion of its effect on i;oo (lower panel). The horizontal axis 
gives the peak- value, gf^'^'^^, of the bump in ghir) in the 
subsonic region (i.e. the bumps in Fig. 1). 

of the bump (Ar c± 0.0015i?*). The line acceleration with 
the extra bumps is shown in Fig. ||. 

The solution of the momentum equation, with the con- 
dition that it passes smoothly through the sonic point, 
gives the velocity at the lower boundary and hence the 
mass-loss rate. The upper panel of Fig. || shows the re- 
sulting mass-loss rates as a function of the peak value of 
the bump in the line acceleration in the subsonic region. 
We see that as the line acceleration in the subsonic region 
increases, M increases. 

2.3. The effect of an increased M on Voo 

Once M is fixed by the processes in the subsonic region, 
the radiative acceleration in the supersonic region then 
determines the terminal velocity v^o that the wind will 
reach. This can easily be seen in the following way. Inte- 
grating the momentum equation (Eq. ^ in the supersonic 
region from the critical point Tc to infinity, and ignoring 
the influence of the gas pressure, gives 



f°° 11 
J 9h{r) dr = - Vesc^ + 2 "° 



so 



2 / gi^{r) dr - Vc 



(6) 



(7) 



Here we have used the observed property that Voo 3> a 
and that Tc — ro <C i?* , so Vc — R*. Eq. |7| says that Voo 
is determined by the integral of gh{r) in the supersonic 
region. 

The radiative acceleration in the supersonic part of the 
wind will decrease as M is forced to increase by an increase 
in the radiative acceleration in the subsonic part of the 



wind. This is because the optical depth of the optically 
thick driving lines, which is proportional to the density 
in the wind, will increase. Thus an increase in M results 
in an increase of the line optical depth. This results in a 
decrease of gh in the supersonic region, which gives a lower 
terminal velocity Voo- We will estimate this effect below. 

Assume that the radiative acceleration by lines de- 
pends on the optical depth in the wind, as given by CAK 
theory (Castor et al. 1975). 



gi^ir) = 5e M(t) 



47^7-2, 



kt~ 



(8) 



where k and a are constants and ge is a reference value 
describing the acceleration due to electron scattering. It 
is given by gc = a'^'L ■ The optical depth parameter is 



t 



cTeVthPidr/dv) 



(9) 



where Vth is the mean thermal velocity of the protons. 
Let us define 5™''(r') as the radiative acceleration in the 
supersonic part of the initial wind model, i.e. without the 
increased mass-loss rate due to the bump in the subsonic 
region, and gh (r) as the radiative acceleration of the model 
with the increased M. From Eqs. |8| and ^ with Eq. || we 
find that 



5r* 



r'^vdv/dr 
[r'^vdv / dr)™ 



M 



,init ' 



M 



(10) 



where the superscript "init" refers to the initial model. 

Let us now compare the terminal velocities of the ini- 
tial model without the bump, to that with the increased 
mass-loss rate due to the bump, in a simple but crude 
way, by solving the momentum equation in the supersonic 
part of the wind. If we neglect the terms due to the gas 
pressure and due to the gravity, the momentum equation 
in the supersonic part of the wind reduces to 



dv 

V— ~ gh{r) 
dr 



(11) 



Solving the equation for the initial model and the model 
with the increased M results in the following expression 



dv 
dr 



dv 
dr 



M 



M 



(12) 



So the ratio between the terminal velocities of the models 
with and without the increased mass-loss rate is 



• init ^1 a/(2-2Q) 
M 



M 



-3/4 



M 



M 



init 



(13) 



where we adopted a — 0.60 (Pauldrach et al. 1986) for 
the last expression. We see that v^o will decrease roughly 

as M '^^^ when the mass-loss rate increases. The result is 
shown in the lower panel of Fig. 0. 
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We realize that this numerical test is a drastic sim- 
plification of the real situation: (a) we have assumed an 
isothermal wind; (b) we have taken the lower boundary 
at a fixed density; (c) we have ignored possible changes in 
the ionization of the wind due to changes in M and (d) 
we have ignored the role of the gas pressure and of grav- 
ity in estimating the change in Voc- However, this simple 
test serves the purpose of explaining qualitatively that the 
mass-loss rate depends on the radiative acceleration in the 
subsonic part of the wind only, and that an increase in the 
mass-loss rate due to an increase of (/l in the subsonic re- 
gion will also be accompanied by a decrease in Vao ■ In the 
rest of the paper, we will quantitatively calculate radiative 
accelerations and mass-loss rates with a method which will 
be described in Sect. |[ 

Thus, an increase in the radiative acceleration in the 
subsonic region of the wind results in an increase of M and 
a decrease in Voo- So, in order to understand the origin 
of the bi-stability jump of radiation driven winds, and 
to predict its effect on M and Voo, we should pay close 
attention to the calculated radiative acceleration in the 
subsonic part of the wind. 

3. The method to predict M 

In order to understand the nature of the bi-stability jump, 
we calculate a series of radiation driven wind models for 
supergiants in the range of Tcff = 12 500 to 40 000 K. The 
calculation of the radiative acceleration of the winds re- 
quires the computation of the contributions of a very large 
number of spectral lines. To this end, we first calculate the 
thermal, density and ionization structure of a wind model 
computed with the non-LTE expanding atmosphere code 
ISA-WIND (de Koter et al. 1993)(for details, see Sect. ^. 
We then calculate the radiative acceleration by following 
the fate of a very large number of photons that are released 
from below the photosphere into the wind, by means of a 
Monte Carlo technique. In this section, we describe the ba- 
sic physical properties of the adopted Monte Carlo (MC) 
technique which was first applied to the study of winds 
of early-type stars by Abbott & Lucy (1985). Then, we 
describe the calculation of the radiative acceleration by 
lines with the MC method, and finally the method for 
calculating theoretical mass-loss rates. 

3.1. Momentum transfer by line scattering 

The lines in the MC method are described in the Sobolev 
approximation. This approximation for the line acceler- 
ation is valid if the physical conditions over a Sobolev 
length do not change significantly, i.e. 



1 
7 



df , 1 
f l«- 
ar Vt 



dv 
dr 



(14) 



where / is any physically relevant variable for the line driv- 
ing, e.g. density, temperature or ionization fraction. Vt is 



a combination of thermal and turbulent velocities. Eq. 
shows that the validity range of the Sobolev approxima- 
tion is in practice somewhat arbitrary, since it depends 
on the value of the turbulent velocity which is poorly 
known. Nevertheless, the Sobolev approximation is often 
used (e.g. Abbott & Lucy 1985) and we will also adopt it 
in calculating the line acceleration and mass loss, mainly 
because of computational limitations. We cannot exclude 
that due to the use of the Sobolev approximation we may 
predict quantitatively inaccurate line accelerations below 
the sonic point. However, if an exact treatment would be 
followed, then this is expected to have a systematic ef- 
fect on the line acceleration for all models. Therefore, we 
do not expect our conclusions regarding the origin of the 
bi-stability jump to be affected. 

The Sobolev approximation implies that for scatter- 
ings in the frame co-moving with the ions in the wind 
(co- moving frame, CMF), the incident and emerging fre- 
quencies are both equal to the rest frequency of the line 
transition i'q in the CMF. 



(15) 



where vl^^ and ly'^^^ are the incident and emerging frequen- 
cies in the CMF. In terms of quantities seen by an outside 
observer, these two CMF frequencies are given by: 



(1 



and 



(16) 



(17) 



where fin and Vont are the incident and emergent frequen- 
cies for an outside observer; /lin and /lout are the direction 
cosines with respect to the radial flow velocity of the pho- 
tons at the scattering point and v is the radial flow velocity 
of the scattering ion for an outside observer. Thermal mo- 
tions of the scattering ions are assumed to be negligible 
compared to the motion of the outward flow. Note that we 
adopted the same velocity v for the ion before and after 
the photon interaction (Eqs. |l^ and [T^). This is justifled 
since the change in velocity due to the transfer of momen- 
tum from a photon to an ion is very small, i.e. about 10^ 
cm s^^ per scattering.. Therefore, the change in and 
Vout is mainly determined by a change in direction angle. 

Combining Eqs. |l^ and |l^ gives the conservation of co- 
moving frequency in a scattering event (Abbott & Lucy 
1985). 



/ , Mill V \ / -I 

t'in ( 1 - ) = I'out ( 1 



A^out V 



(18) 



Because the energy and momentum of a photon are E — 
hu and p — hv j c, the equation can be rewritten in the 
following way: 



Pin Mill Pout Mout 



(19) 
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Eq. links the change in radial momentum of a photon 
in an interaction with an ion with velocity v to the energy 
loss of the photon. In order to determine the line accel- 
eration gi^ we will need to derive the momentum transfer 
from the photons to the ions in the wind. 

For an outside observer, the conservation of radial mo- 
mentum is; 



mvi 



Mil 



mv2 



Mout 



(20) 



where m is the mass of the moving ion and vi and V2 are 
the radial velocities of the ion just before and after the 
scattering. For an outside observer: 



I^in = I'D (1 + 

and 

Vout = (1 + 



Min^'^ 



(21) 



(22) 



Again, the change in frequency is dominated by the change 
in direction angle. So the change in radial velocity per 
scattering, Aw ~ W2 ~ ^^i, is small compared to v and is 
given by 



Aw = 



V2 - Vi 



Mil 



Since w ^ c, Eq. 



Av 



V2 - Vi 



-)Min - 
c 

\ becomes 

(Min - 



(1 



MoutV 



Moutj 



)Mout(23) 



(24) 



This relation describes the velocity increase of the ion de- 
pending on the directions /ii„ and /iout of the photon. In 
case /Zin = Mout then Aw = 0, as one would expect. The 
increase in the radial momentum Ap = m Aw of the 
scattering ion is now given by: 



Ap = m{v2 — vi) 



c 

hv\r 



(Mil 



Moutj 



AE 

V 



(25) 



where AE — Ein — i?out is the loss of radiative energy. 
This equation shows that the increase in the momentum 
of the ions can be calculated from the loss of energy of the 
photons when these are followed in their path through the 
wind by means of the Monte Carlo method. 

Multiplying both sides of Eq. 25 by v and using the 
fact that for each scattering W2 — wi so w ~ [vi + W2)/2, 
gives: 



1 



m{V' 



(26) 



Equation [2^ says that the gain of kinetic energy of the 
ions in the radial direction equals the energy loss of the 
photons. 



3.2. Scattering and absorption processes in the MC 
calculations 

The radiative acceleration as a function of distance is cal- 
culated by means of the MC technique by following the 
fate of the photons using the program MC-WIND (de Koter 
et al. 1997). In the calculation of the path of the photons 
we have properly taken into account the possibility that 
the photons can be scattered or absorbed & re-emitted 
due to true absorption or eliminated because they are scat- 
tered back into the star. 

The radiative transfer in MC-wind is calculated in the 
Sobolev approximation. Multiple line and continuum pro- 
cesses are included in the code. The continuum processes 
included are electron scattering and thermal absorption 
and emission. The line processes included are photon scat- 
tering and photon destruction by coUisional de-excitation. 
In deciding whether a continuum or a line event takes 
place, we have improved the code in the following way: 
The key point of the Monte-Carlo "game" is that line pro- 
cesses can only occur at specific points in each shell of the 
stellar wind, whereas continuum processes can occur at 
any point. The correct way of treating the line and con- 
tinuum processes is by comparing a random optical depth 
value to the combined optical depth for line and contin- 
uum processes along the photon's path. First, this com- 
bined optical depth is compared to a random number to 
decide whether a continuum or a line process takes place. 
This first part of the treatment is basically the same as de- 
scribed by Mazzali & Lucy (1993) for the case of line and 
electron scatterings only. Now, after it has been decided 
that the process will be a continuum process, a second ran- 
dom number is drawn to decide which continuum process 
will take place, an electron scattering or absorption. 

3.3. The calculation of the radiative acceleration gi^{r) 

The radiative acceleration of the wind was calculated by 
following the fate of the photons emitted from below the 
photosphere with the MC technique. To this purpose the 
atmosphere is divided into a large number of concentric, 
thin shells with radius r, thickness Ar containing a mass 
Am(r). 

The loss of photon energy due to all scatterings that 
occur within each shell are calculated to retrieve the total 
line acceleration g-L^r) per shell. The total line acceleration 
per shell summed over all line scatterings in that shell 
equals 



1 S Ap{r) 
Am(r) Ai 



(27) 



where p{r) is the momentum of the ions in the shell. The 
momentum gained by the ions in the shell is equal to the 
momentum lost by the photons due to interactions in that 
shell. Using the relationship between Am{r) and Ar for 



Jorick S. Vink et al.: The nature of the Bi-stability jump 



7 



thin concentric shells, Am(r) = 47rr^p(r)Ar, and the de- 
rived relation between momentum and energy transfer of 
the photons Ap — AE/v (Eq. g^ir) can be rewritten 
as 



5lW 



1 



S AE{r) 



Anr^ p{r)Ar v{r)At 



(28) 



where Y,AE{r) is sum of the energy loss of all the pho- 
tons that are scattered in the shell. Now using mass con- 
tinuity (Eq. H) and the fact that the total energy trans- 
fer E AE{r) divided by the time interval At equals the 
rate at which the radiation field loses energy, —AL{r), i.e. 
E AE{r)/At ~ — AL(r), the expression for ghir), which 
is valid for each shell, simply becomes (Abbott & Lucy 
1985) 



1 AL{r) 
M At 



(29) 



The line list that is used for the MC calculations consists 
of over 10^ of the strongest lines of the elements H - Zn 
from a line list constructed by Kurucz(1988). Lines in the 
wavelength region between 50 and 7000 A are included 
in the calculations with ionization stages up to stage vi. 
Typically about 2 10^ photon packets, distributed over the 
spectrum at the lower boundary of the atmosphere were 
followed for each model, i.e. for each adopted set of stellar 
and wind parameters. For several more detailed models 
we calculated the fate of 2 10^ photon packets. The wind 
was divided in about 50-60 concentric shells, with many 
narrow shells in the subsonic region and wider shells in 
supersonic layers. The division in shells is essentially made 
on the basis of a Rosseland optical depth scale. Typical 
changes in the logarithm of this optical depth are about 

0. 13. 

3-4-. The determination of AI 

We predict the mass-loss rates for a grid of model atmo- 
spheres to study the behaviour of M near the bi-stability 
jump. For a given set of stellar parameters we calculate 
the mass loss in the following way: 

1. For fixed values of L, T^tf, R* and M^s we adopt several 

values of the input mass loss M ^ (within reasonable 
bounds predicted by CAK theory). 

2. For each model we adopt a wind with a terminal veloc- 
ity of 1.3, 2.0 or 2.6 times the effective escape velocity, 
given by 



2GMeff 



(30) 



A P-type velocity law with f3 — I was adopted, ap- 
propriate for OB stars (Groenewegen & Earners 1989; 
Puis et al. 1996) 



3. For each set of stellar and wind parameters we calcu- 
late a model atmosphere with ISA- wind (see Sect. ||). 
This code gives the thermal structure, the ionization 
and excitation structure and the population of the en- 
ergy levels of all relevant ions. 

4. For each model the radiative acceleration was calcu- 
lated with the MC-wiND program that uses the Monte 
Carlo method described above. 

5. For each set of stellar parameters and for each adopted 
value of we check which one of the adopted mass- 
loss rates is consistent with the radiative acceleration. 
This consistency was checked in the following way: 
Neglecting the term due to the gas pressure, one can 
write the equation of motion in the following way: 



dv 
dr 



GM.ff 



(31) 



Using the expression for the line acceleration (Eq. p9[ ) 
and integrating the equation of motion (Eq. ^l[ ) from 
the stellar surface to infinity gives 



2 M (woo' 



= AL = M gL{r)dr (32) 
Jr, 



AL — EAi(r), is the total amount of radiative energy, 
summed over all the shells, that is lost in the process 
of line-interaction and is transfered into kinetic energy 
of the ions as given in Eq. Equation ^ states that 
the momentum transfered from the radiation into the 
wind is used to lift the mass loss out of the potential 
well and to accelerate the wind to v^o ■ Only one value 
of M will satisfy this equation (Lucy & Abbott 1993). 
This is the predicted mass-loss rate. 

We note that Eq. ^ only describes the "global" consis- 
tency of the mass-loss rate with the radiative acceleration. 
For the set-up of the model atmosphere the velocity law 
v{r) is needed as input. This means that although the M 
calculation is globally consistent in terms of kinetic wind 
energy, the velocity is not necessarily locally consistent, 
since the equation of motion is not solved. Instead, we 
have used observed values for Voo and /3 for the velocity 
law. Since the total amount of radiative energy in Eq. |3j 
is mainly determined in the supersonic region, where the 
Sobolev approximation is an excellent approximation, AL 
is accurately calculated. This implies that if one adopts 
the correct values for the terminal velocity, one may pre- 
dict accurate values for Ml 

4. The model atmospheres 

The calculation of the mass-loss rates by the method de- 
scribed in the previous section requires the input of a 
model atmosphere, before the radiative acceleration and 
M can be calculated. 

The model atmospheres used for this study are calcu- 
lated with the most recent version of the non-LTE unified 
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Improved Sobolev Approximation code ISA-wind for stars 
with extended atmospheres. For a detailed description of 
this code we refer to de Koter et al. (1993, 1997). Here, 
we just make a few relevant remarks. 

ISA-WIND treats the atmosphere in a unified manner, 
i.e. no artificial separation between photosphere and wind 
is assumed. This is distinct from the so-called "core-halo" 
approaches. In the photosphere the density structure fol- 
lows from a solution of the momentum equation taking 
into account gas and radiative pressure on electrons. The 
velocity law follows from this density structure via the 
mass continuity equation. Near the sonic point, a smooth 
transition is made to a /3-type velocity law for the super- 
sonic part of the wind (see Eq. ||). 

The temperature structure in the wind is computed 
under the assumption of radiative equilibrium in an ex- 
tended grey LTE atmosphere. The temperature in the 
wind is not allowed to drop below a certain minimum value 
Tmin = 1/2 Tcff (Drew 1989). Finally, the chemical species 
included explicitly in the non-LTE calculations are H, He, 
C, N, O and Si. The complexity of the model atoms is 
similar to that used by de Koter et al. (1997). For the iron- 
group elements, which are important for the radiative ac- 
celeration, we calculate the ionization/excitation equilib- 
rium in the modified nebular approximation (see Schmutz 
1991). In this representation the ionization equilibrium is 
given by 



{{1-Ow+ow (^^) 



1/2 



LTE 



(33) 



where and To are the electron density and tempera- 
ture, Nj and iVj+i are the ion population numbers, Tr 
= TR(r, j) is the radiation temperature of ion j at ra- 
dial depth r, and is a geometrical dilution factor as 
defined by Schmutz et al. (1990). The last factor of Eq. 
|33| is the LTE ionization ratio for a temperature T^(r,j). 
The parameter ^, introduced by Abbott & Lucy (1985), 
represents the fraction of recombinations going directly 
to the ground state. The values of TR(r, j) are obtained 
by inverting the above equation, using all 19 ionization 
ratios available from the ISA- wind calculation. The radi- 
ation temperature of an explicit ion is used that has its 
ionization potential closest (but lower) to that of the metal 
ion of interest. For instance, the N ii/iii ratio is used to 
define the ionization equilibrium of Fe iii/iv. 

The excitation state of metastable levels is assumed to 
be in LTE relative to the ground state. For all other levels 
we adopt "diluted" LTE populations, defined by 



ni 



W 



Hi 



LTE 



(34) 



where riu and ni are the excitation population numbers for 
the upper and lower levels. Clearly, the simplified treat- 
ment of the iron-group metals is prompted by the com- 
putationally intensive nature of the problem at hand. It 



needs to be improved in the future, but we do not ex- 
pect that our conclusions regarding the nature of the bi- 
stability jump would be affected. (We return to this in the 
discussion in Sect. ||). 

5. The predicted bi-stability jump 



Using the procedure as described in Sect. 3.4, we cal- 



culated mass-loss rates for stars with a luminosity of 
L* = 10^ Lq and a mass of Af* = 20Mq. The models 
have effective temperatures between 12 500 and 40 000 K 
with a stepsize of 2500 K. These parameters are approx- 
imately those of OB supergiants, for which Lamers et al. 
(1995) found the bi-stability in v^. 

We calculated M for wind models with a /3-type ve- 
locity law with (3 = 1 (Eq. ^ for three values of the 
Voo/vcsc = 2.6, 2.0 and 1.3. Lamers et al. (1995) found 
that Voo/vcsc — 2.6 for stars of types earlier than Bl, 
and Voo/vcsc — 1-3 for stars of types later than B2. 
For the determination of Vcsc we used the effective mass 
Mcs = 17.4 Mq, with Te ^ 0.130. 

The stellar parameters for the calculated grid are in- 
dicated in Table The models are calculated for solar 
metallicities. 

5.1. The predicted bi-stability jump in M 

The results are listed in Table 0. This Table gives the 
values of Tcff, i?*, Ucsc and M for each temperature and 
for the three values of Voo/vcsc- We also give the value of 
the wind efficiency factor 77, which describes the fraction 
of the momentum of the radiation that is transferred to 
the ions 



Mvo 



(35) 



The fraction of the photon energy that is transferred into 
kinetic energy of the ions is also listed (in column 8). The 
values for this energy efficiency number l^L/ L are a factor 
of about 10"'^ smaller than the wind momentum efficiency 
number rj, which is given in column (7). This is because a 
photon transfers a large fraction of its momentum during a 
scattering, but only a very small fraction (of order v/c) of 
its energy. The last column of Table |l| marks three models 
that will be discussed in more detail in Sect. ||. 

The results are plotted in Fig. ||. For each of the three 
values of Voo/vesc the value of M is decreasing for decreas- 
ing Tcff between 40 000 and 30 000 K and also between 
22 500 and 12 500 K. Between about Tcff = 27 500 K and 
Tcff = 20 000 K (slightly dependent on Voo/vesc) the mass 
loss increases with decreasing Teff. These increments in 
M roughly coincide in Tcff with the observed bi-stability 
jump in Voo/vesc near spectral type Bl, at about 21 000 
K. For the ratio of Vqo/vosc = 2.6, the increase in M be- 
tween model A and B equals 45 %. We know from the 
observations that Vqo/vcsc jumps from 2.6 at the hot side 
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Table 1. Stellar parameters of the grid of calculated models. 

log L/Lq = 5.0, M ^ 2OM0, Te = 0.130, M^s = 17.4 M©, /3 = 1, solar mctallicity. 



Voo/vosc Tcff R* Vesc Voo log M 7] AL / L model 



(K) (Rq) (kms-^) (kms-^) (Mp/yr) (in 10"^) 



12 500 


67.7 


310 


410 


- 6.32 


0.095 


0.103 


15 000 


47.0 


380 


490 


- 6.39 


0.097 


0.126 


17 500 


34.5 


440 


570 


- 6.28 


0.146 


0.221 


20 000 


26.4 


500 


650 


- 6.22 


0.192 


0.332 


22 500 


20.9 


560 


730 


- 6.15 


0.254 


0.493 


25 000 


16.9 


630 


810 


- 6.12 


0.302 


0.653 


27 500 


14.0 


690 


900 


- 6.40 


0.174 


0.414 


30 000 


11.8 


750 


980 


- 6.58 


0.126 


0.326 


32 500 


10.0 


810 


1060 


- 6.58 


0.136 


0.382 


35 000 


8.6 


880 


1140 


- 6.43 


0.207 


0.626 


37 500 


7.6 


940 


1220 


- 6.37 


0.255 


0.826 


/in nnn 


D.D 


1 nnn 
iUUU 


1 Qnn 
ioUU 


- D.ZD 


U.ooU 


i.ziu 


12 500 


67.7 


310 


630 


- 6.74 


0.056 


0.073 


15 000 


47.0 


380 


750 


- 6.62 


0.088 


0.138 


17 500 


34.5 


440 


880 


- 6.49 


0.139 


0.254 


20 000 


26.4 


500 


1000 


- 6.41 


0.191 


0.398 


22 500 


20.9 


560 


1130 


- 6.32 


0.264 


0.620 


25 000 


16.9 


630 


1250 


- 6.48 


0.203 


0.530 


27 500 


14.0 


690 


1380 


- 6.73 


0.125 


0.360 


30 000 


11.8 


750 


1500 


- 6.76 


0.128 


0.400 


32 500 


10.0 


810 


1630 


- 6.71 


0.155 


0.527 


35 000 


8.6 


880 


1750 


- 6.59 


0.220 


0.801 


37 500 


7.6 


940 


1880 


- 6.57 


0.247 


0.969 


40 000 


6.6 


1000 


2000 


- 6.48 


0.325 


1.356 


12 500 


67.7 


310 


810 


- 6.95 


0.045 


0.070 


15 000 


47.0 


380 


980 


- 6.85 


0.067 


0.126 


17 500 


34.5 


440 


1140 


- 6.69 


0.114 


0.248 


20 000 


26.4 


500 


1300 


- 6.54 


0.184 


0.458 


22 500 


20.9 


560 


1460 


- 6.59 


0.184 


0.517 


25 000 


16.9 


630 


1630 


- 6.79 


0.129 


0.403 


27 500 


14.0 


690 


1790 


- 6.95 


0.098 


0.337 


30 000 


11.8 


750 


1950 


- 6.92 


0.115 


0.430 


32 500 


10.0 


810 


2120 


- 6.86 


0.143 


0.579 


35 000 


8.6 


880 


2280 


- 6.76 


0.194 


0.845 


37 500 


7.6 


940 


2440 


- 6.71 


0.233 


1.089 


40 000 


6.6 


1000 


2600 


- 6.68 


0.266 


1.327 



c 



of 21 000 K to 1.3 at the cool side of 21 000 K (Lamers et 
al. 1995). Including this observed jump in Voo/vesc in the 
mass-loss predictions, provides an even steeper increase in 
M from models A and B to the smaller value of Voo/vesc 
= 1.3, as is shown in the lower part of Fig. This figure 
shows an increase in M of about a factor of five between 
Teff = 27 500 and 20 000 K. This is our prediction for a 
bi-stability jump in M. 

The exact position of Teff of the bi-stability jump in 
Fig. ^ is somewhat ambiguous, since Woo is adopted from 
observations, and does not directly follow from our mod- 
els. For a discussion on the exact position of the jump in 
Teff, see Sect. ||. 

To test the sensitivity of our predictions of mass-loss 
rates for different shapes of the velocity law, we calculated 



another series of models with /3 = 1.5 . Since the differ- 
ences are only about 10 %, we conclude that the predicted 
mass-loss rates are only marginally sensitive to the shape 
of the adopted velocity law. 

5.2. The predicted bi-stability jump in rj 

Another view at these results can be obtained by plotting 
the wind efhciency factor rj. Figure^ shows the behaviour 
of 77 as a function of Teff for the same grid of models as 
was presented for the mass-loss rates in the upper panel 
of Fig. I 

Fig. ^ clearly shows that rj is not a constant func- 
tion of Toff- The overall picture shows that for the three 
values of Voo/vosc, V decreases as Toff decreases. This is 
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Fig. 3. Upper panel: The calculated mass-loss rates M as a function of Tgs for three values of the ratio Voo/vesc- The 
values for Vqo/vosc are indicated in the lower left corner. The stellar parameters arc log L/Lq = 5.0, M = 20 Mq 
and P = 1.0; all models are calculated for solar metallicities. Lower panel: The predicted bi-stability jump in M from 
models with the observed ratios of Voo/vesc = 2.6 for T^s > 21,000 K and Vcxa/vesc = 1-3 for T^g < 21,000 K, as 
indicated in the lower left corner. 



probably due to the fact that the maximum of the flux 
distribution shifts to longer wavelengths. At A > 1800 
A there arc significantly less lines than at A < 1800 A. 
Therefore, radiative acceleration becomes less effective at 
lower Teff. In the ranges of 40 000 < T^s < 30 000 and 
20 000 < Teff < 12 500 K, r] is almost independent of the 
adopted value for Voo/^esc- This means that the behaviour 



of 1] is intrinsically present in the model calculations and 
does not depend on the values adopted for Voo/vesc- 

In the range of 30 000 < Tos < 20 000 K, the situation 
is reversed, rj now increases by a factor of 2 to 3. This 
means that the wind momentum loss, Mvoo is not constant 
over the jump, but instead, jumps by a factor of 2 - 3 also. 
Since Voo drops by a factor of about two, M is expected to 
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Fig. 4. The wind efficiency number i] = Mvoo/ {L^,/c) as a function of Tcs for three values of the ratio Voo/vcsc- These 
values are indicated in the upper right corner. Note the steady decrease of 77 to lower temperatures, except the jump 
of about a factor 2 or 3 near 25 000 K. 



jump by a factor of about five, which was already shown 
in the lower panel of Fig. ^. 

The behaviour of 77 as a function of Teg is not exactly 
the same for the three different series of models. First, 
the size of the jump is different. Second, the jump occurs 
at somewhat different temperatures. This is no surprise, 
since the ionization equilibrium does not only depend on 
T, but on p as well, a smaller value of the velocity Vqo, 
means a larger density p in the wind. Hence, the jump 
is expected to start at a larger value of Tgff for a smaller 
value of Voo/vesc- This behaviour for the position of Tog 
of the jump can be seen in M in Fig. |^ and in rj in Fig. ^. 

6. The origin of the bi-stability jump 

In the previous section we have shown that the mass-loss 
rate increases around Tee = 25 000 K. The next step is to 
investigate the physical process that causes the bi-stability 
jump. Therefore, we will look into the details of the line 
acceleration gi^ (r) for three models around the bi-stability 
jump. For these models (A, B and C in Table |l| and Fig. ||) 
we made improved Monte-Carlo calculations, using 2 x 
10^ packets of photons, to derive more details about the 
radiative acceleration. 

First, we will investigate the line acceleration ghir) of 
the model at the hot side of the bi-stability jump. This 
model A with T^ff = 27 500 K and Voo/vesc = 2.6, is our 
basic model. Then, we will compare model A to model B 
that has the same V(x>/vcsc, but is situated on the cool side 
of the bistability jump, where Tcs = 25 000 K. By com- 
paring models A and B, we can investigate the intrinsic 



increase in M of 45 % in our model calculations due to 
the lower Tcb- The next step is to compare ghif) of model 
B and model C which also has Tcs — 25 000 K, but a 
smaller ratio Voo/vesc = 1-3. By comparing model B and 
C, we can obtain information about the effects of a jump 
in Vac- Finally, we check our approach for self-consistency 
by simultaneously calculating the mass-loss rate and ter- 
minal velocity. 

6.1. The main contributors to the line acceleration 

Model A has a mass- loss rate of log M — -6.95. The be- 
haviour of the line acceleration as a function of the dis- 
tance from the stellar surface, ghir) is shown in Fig. ^. The 
sonic point is reached at a distance of 1.025 i?,. It is clear 
that most of the line driving is produced far beyond the 
sonic point. But, as was explained in Sect. the impor- 
tant region that determines the mass-loss rate is below the 
sonic point. Therefore, the part of the atmosphere around 
the sonic point is enlarged in Fig. ^(b). 

To investigate the origin of the jump, it is useful to 
know which elements are effective line drivers in which 
part of the stellar wind. Therefore, extra Monte-Carlo 
calculations were performed. The first extra Monte-Carlo 
simulation was performed with a line list containing only 
Fe lines. The second one was performed with a line list 
containing the lines of the elements C, N and O. 

Figure ^(b) shows that Fe is the main line driver below 
the sonic point. C, N and O, are important line drivers in 
the supersonic part of the wind, which can be seen in ^(a). 




Fig. 5. The line acceleration of model A {Tcft— 27 500 K and Voo/vosc = 2.6), from 1 to 15 i?* (left) and around 
the sonic point (right), (a) The solid line shows the total gi^ as a function of the distance. The dashed line is the 
contribution by C, N and O only. The dotted line shows the contribution by Fe lines. Some values for the velocity are 
indicated on top of the figure, (b) The region around the sonic point is enlarged. The sonic point is reached at x = 
1.025. Note the bump in the gh{r) just below the sonic point, which is largely due to Fc lines. 



C, N and O contribute roughly 50 % of the line acceler- 
ation in the supersonic part of the wind. Not indicated 
here, but relevant to mention is that Si, CI, P and S are 
other important line drivers in the supersonic part of the 
wind. Ni was found not to be an important line driver in 
any part of the stellar wind at all. 

The mass-loss rate is determined by the radiative ac- 
celeration below the sonic point, and the terminal velocity 
is determined by the radiative acceleration in the super- 
sonic part of the wind. So our results show that the mass- 
loss rates of hot star winds are mainly determined by the 
radiation pressure due to Fe! The terminal velocities are 
mainly determined by the contributions of C, N and O. 

6. 2. The effect of the Fe ionization 

To understand the origin of the bi-stability jump in M, we 
investigate the line acceleration due to Fe. The ionization 
balance of Fe for models A and B is plotted in Fig. |^, top 
and bottom respectively. The right hand figures show the 
enlargement of the ionization balance in the region near 



the sonic point. In Model A {T^s=27 500 K) Fe V has a 
maximum around x = 1.004, which can be seen in Fig. |^ 
(b). Then, due to the outward decreasing temperature, 
Fe V decreases in favour of Fe IV, which peaks around 
X — 1.008. Next, one may expect Fe iv to decrease in 
favour of Fe ill. However, around x — 1.013 Fe IV re- 
ionizes due to a decrease of the density p. In this region 
of the atmosphere, where dv/dr is rapidly increasing, the 
effect of the decreasing p is larger than the effect of the 
decreasing T. 

Fig. H (b) clearly shows that Fe IV is the dominant ion- 
ization stage in the subsonic region of the stellar wind. In 
the region just below the sonic point, the ionization frac- 
tion of Fe IV is 90 - 100 % whereas that of Fe ill is less than 
10 %. However, this does not necessarily mean that Fe iv 
is the main line driver. To investigate the contribution to 
the line acceleration of the different ionization stages of 
Fe some extra Monte-Carlo simulations were performed. 
One simulation included only the lines of Fe ill, another 
simulation included just the lines of Fe iv. The results for 
(7l for Fe ill and Fe iv are plotted in Fig. 0. 
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Fig. 6. The ionization fraction of Fe as a function of distance. The upper panels are for model A and the lower panels 
for model B. (a) Fe ionization for model A from x = 1 to 15. (b) Model A, enlarged around sonic point, (c) Fe 
ionization for model B from x = I to 15. (d) Model B, enlarged around sonic point. 



It is surprising to note that, although Fe iv is the dom- 
inant ionization stage throughout the wind, most of the 
driving is contributed by Fe ill. Below the sonic point 
Fe III is clearly the most important iron line driver (see 
Fig. 0(b)). 

From the data shown in Figures ^ and ^ we con- 
clude that the mass-loss rate of winds from stars with 
Tcif ~ 27 500 K is mainly determined by the radiative 
acceleration due to Fe ill lines. This suggests that the bi- 
stability jump is mainly due to changes in the ionization 
balance of Fe. We test this hypothesis in the next section. 

6.3. The effect ofT^s on M 

In the previous section we have shown that the mass loss 
of model A is dominated by radiative acceleration due to 
Fe III lines. In this section we investigate changes in the 
radiative acceleration due to Fe as Tcs decreases. This may 
explain the increase of M near the bi-stability jump. To 
this purpose we compare the ionization and (?l of models 
A and B in detail. 

The ionization balance of model B is shown in Fig. |^(c) 
and (d). It shows that, due to a lower temperature, the 



decrease of the Fe iv fraction drops to smaller values than 
for model A, which was shown in Fig. ^(b). The ionization 
fraction of Fe ill below the sonic point in the case of model 
B is up to almost 40 %. To see whether this extra amount 
of Fe III can cause the increase in the line acceleration, we 
must look at of Fe for model B. 

Since model A and B have different Tcs at the same , 
they have a different radiative surface flux. The radiative 
acceleration will be proportional to this flux. In order to 
compare the values of of the two models, we scale the 
results to a flux of a Tcff = 25 000 K model. So 

»r™ = (^)' (36) 

Since Tcff^ for constant luminosity, this is also a 

scaling to the Newtonian gravity of the models. 

Figure ^ shows the normalized gh of Fe for the models 
A (top) and B (bottom). The right hand figures show an 
enlargement of the region near the sonic point. It shows 
that for model B gi^ of Fe ill around the sonic point is more 
than a factor two larger than for model A (see Figs, ^(b) 
and (d)). This extra amount of Fe ill in model B causes 
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Fig. 7. The contribution of several Fc ions to gi, as a function of distance from the stellar surface for model A. (a) The 
full distance range of a; = 1 to 15. (b) The region around the sonic point is enlarged. The legend indicates the ionization 
stage. Some values for the velocity are indicated on the top of the figure. Note that the strongest contribution to gL 
below the sonic point is due to Fe ill, although the ionization fraction of this ion is less than 10 %. 



an increase in the total in the subsonic part of the wind 
also, as can be seen in Fig. |9|(b). 

We conclude that the increase in mass loss from model 
A to B is due to the larger radiative acceleration ( compared 
to the gravity) of model B by a larger ionization fraction 
of Fe III below the sonic point. 

6.4- The effect of Voa 

Now the effect of on v^a will be examined. Therefore, 
Model B is compared to model C. We remind that models 
B and C have the same T^s, and hence the same radiative 
flux and gravity, but model C has a twice as small value of 
Voo/vesc as model B. Figure ^(a) shows the normalized (?l 
for models A, B and C. As expected, gh{r) for model C is 
significantly smaller than gi,{r) for models A and B. This 
is obviously due to the smaller value of Voo- The integral 
J dhir) dr in Fig. ||(a) for model A and B is larger than 
for model C. The values of / ghir) dr for the models are 
2.34 X 10^^ and 1.92 x lO^^ cm^ s'^ for models A and B 
respectively, and 6.12 x 10^^ cm^ s^^ for model C. Using 
Eq. and the values of Wesc from column (4) in Table |l|, 
the output values for v^o can be obtained from the val- 



ues of the integral of gi^. The derived output values for 
Voo for the models are Voo = 2050, 1860 and 920 km 
respectively for the models A,B and C. These values are 
equal within 10 % to the input values for Voo which were 
indicated in column (5) of Table |l|. We can conclude that 
a smaller value for Voc is indeed consistent with a smaller 
value of the integral / gh{i") dr. However, this is not an 
independent check, since the calculated line acceleration 
of optically thick lines (in the Sobolev approximation) is 
inversely proportional to the Sobolev optical depth which 
is proportional to {dv/dr)~^. Hence, assuming a smaller 
terminal velocity will automatically result in a smaller cal- 
culated line acceleration. 

6.5. A self-consistent solution of the momentum equation 

In earlier sections we have demonstrated that the mass 
loss around the bi-stability jump increases. As we have 
used observed values for the ratio Voo/vesc in our model 
calculations, we have not yet provided a self-consistent 
explanation of the observed bi-stability jump in Voo/vesc- 
As a consistency test of our calculations and an attempt 
to explain the observed jump in the ratio Voo/vcsc, we 
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Fig. 8. Normalized of Fc as a function of distance from the stellar surface for the models A and B. (a) Normalized 
gi^ for the different Fe ionization stages of model A. The legend indicates the ionization stage. Some values for the 
velocity are indicated on the top of the figure, (b) model A, enlarged around the sonic point, (c) Normalized for 
the different Fe ionization stages of model B. (d) model B, enlarged around sonic point. 

Table 2. Force multipliers and consistent models. 



log L/Lq 


= 5.0, M = 


20Mq, Te = 


0.130, Moff 


= 17.4 Mq, /? = 1, solar metallicity. 






Tcs 


(Voo/Vcsc)o 


(Voo/Vosc)l 


(Voo/Vcsc)2 


(Voo/Vcsc)3 CC^'^ k^'^ (Voo/Vosc)4 


log A'/cAK 


log Mmc 


(K) 










{Me/yr) 


(Me/yr) 


17 500 


2.0 


1.5 


1.3 


1.2 0.58 0.2065 1.2 


-6.21 


-6.27 


30 000 


2.0 


2.5 


2.4 


2.4 0.85 0.0076 2.4 


-6.86 


-6.90 



proceeded to solve the momentum equation of line driven 
wind models around the bi-stability jump. The approach 
we take is to combine predicted force multiplier parame- 
ters k and a (see below) from the Monte Carlo calcula- 
tion with the analytical solution of line driven winds from 
CAK. 

We calculated the line acceleration gi^ for several mod- 
els with different Tes using the Monte Carlo method. The 
values of gh were expressed in terms of the force multiplier 
M{t) (Eq. |). Following CAK we tried to express M{t) in 
terms of a power-law fit of the optical depth parameter 
t (Eq. 1). We found that in the range 20 000 < T^s < 
27 500, M{t) is not accurately fit by a power-law, since 
the ionization changes over this critical range in T^g. For- 



tunately, just outside this temperature region, M{t) can 
be accurately represented in terms of k and a, i.e. 

Therefore, we have calculated models with effective tem- 
peratures just below (Toff = 17 500 K) and just above 
(Tcff = 30 000 K) this critical temperature range. Self- 
consistent values of Uoo and M were thus found in the 
following way: 

1. We started with an assumed ratio of Voo/vosc = 2.0 
(See column (2) in Table 0). 

2. The force multipliers M'^(t) were calculated and a 
power-law fit of the type Eq. ^ was derived. The fit 




Fig. 9. (a) The normalized total gh for the models A, B and C as a function of distance. Notice the much smaller 
radiative acceleration in the supersonic region of model C compared to models A and B. (b) An enlargement of 
the region around the sonic point. The sonic point is located around x = 1.025 r/R^,. Notice also the much smaller 
radiative acceleration in the subsonic region of model C compared to models A and B. This is due to the smaller value 
of Vqo/vosc for model C. 



was found to be excellent in the important part of the 
wind between the sonic point and v ~ O.Sdqo- This 
yielded values of a^'^ and k^'^^ . Next, the mass loss 
and terminal velocity were simultaneously calculated 
from these a^'^ and k^'~' parameters using the CAK 
solution of the momentum equation. Note that the so- 
lution with the finite disk correction (Pauldrach et al. 
1986) was not applied, since this is already properly 
taken into account in the values of a'^*-' and k^'^ cal- 
culated in the Monte Carlo technique (see Sect. 
The superscript, MC, to the force multiplier parame- 
ters was added to avoid confusion with k and a for a 
point-like source used by e.g. Kudritzki et al. (1989). 
The ratio Voo/vesc can be derived from the simple CAK 
formulation: 



^esc 



vMC 



I - a 



MC 



(38) 



The value for a^^ for the model of 30,000 K is signif- 
icantly higher than values for a that were calculated 



before (e.g. Pauldrach et al. 1986), since the finite disk 
is already included in the a*^*-^-parameter! 

3. The new calculated terminal velocity ratio Vqo/vcsc 
(column (3) of Table ||) was used in the next iteration. 

4. New mass-loss rates were calculated from the MC ap- 
proach using the procedure as explained in Sect. ^.4[ . 
The mass-loss rates are equal within 15 % to the mass- 
loss rates that can be calculated using the expression 
for M of CAK using a^'^ and fc^^ . 

5. The above procedure (step 1. through 4.) was repeated 
until convergence was reached. After four iterations, 
the ratio Voo/vesc did not change anymore. The inter- 
mediate values of Voo/vcsc are given in columns (3), (4) 
and (5) of Table |^. The final value for the ratio Vqo/vcsc 
is given in column (8). For the hot model (Toff = 30 
000 K) the final ratio Voo/vosc equals 2.4; for the cool 
model (Tcff = 17 500 K) Voo/vosc = 1-2. These values 
are within 10 % of the observed values of Voo/vcsc, i-e. 
2.6 and 1.3 respectively. 

6. CAK mass-loss rates were also calculated from the re- 
sulting final force multiplier parameters fc^^ and 
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(given in columns (6) and (7) of Table || and the final 
mass-loss rates are given in column (9) of this Table. 
Note that the values of M are only marginally differ- 
ent from the mass-loss rates that were calculated from 
the Monte Carlo approach (column (10) of Table 

In summary; we have self-consistently calculated val- 
ues of Voo and M of two models located at either side of 
the bi-stability jump. We have found a jump in terminal 
velocity Voo/vesc of a factor of two, similar as observed by 
Lamers et al. (1995). Moreover, the mass-loss rates calcu- 
lated from the CAK formulation are consistent with those 
obtained from our Monte Carlo approach. This implies 
that the origin of the observed change in the ratio Voo/fcsc 
of a factor of two around spectral type Bl is identical to 
the predicted jump in mass-loss rate of a factor of five due 
to the recombination of Fe IV to Fe III. 

6. 6. Conclusion about the origin of the bi-stability jump 

From the results and figures presented above we conclude 
that the mass-loss rate of early-B supergiants near the bi- 
stability jump is mainly determined by the radiative accel- 
eration by iron. Although Fe iv is the dominant ionization 
stage in the atmosphere of stars near 25 000 K, it is Fe ill 
that gives the largest contribution to the subsonic line ac- 
celeration. This is due to the number of effective scattering 
lines and their distribution in wavelengths, compared to 
the energy distribution from the photosphere. This im- 
plies that the mass-loss rates of B-supergiants are very 
sensitive to the ionization equilibrium of Fe in the upper 
photosphere. Our models show that the ionization frac- 
tion of Fe III increases drastically between Tca= 27 500 
and 25 000 K. This causes an increase in the line accelera- 
tion below the sonic point and in turn increases the mass 
loss near the bi-stability jump. 

7. Bi-stability and the variability of LBV stars 

Luminous Blue Variables (Conti 1984) are massive stars 
undergoing a brief, but important stage of evolution. Dur- 
ing this period they suffer severe mass loss with M values 
of up to lO~^M0yr~-'^. LBVs are characterized by typi- 
cal variations in the order of /SV of 1 to 2 magnitudes. 
Nevertheless, the total bolometric luminosity of the star 
L* seems to be about constant. The reason for the typical 
LBV variations is still unknown. For reviews see Nota & 
Lamers (1997). 

Leitherer et al. (1989) and de Koter et al. (1996) have 
shown that it must be the actual radius of the star that 
increases during these typical variations. Therefore, Tcs 
decreases during the variations, if is about constant. 
In this paper, we have calculated the mass-loss behaviour 
for normal OB supergiants as a function of Teff. Despite 
many differences between OB supergiants and LBVs, we 
can retrieve valuable information about the behaviour of 



M during a typical LBV variation by investigating the 
M behaviour of normal OB supergiants, since both types 
of stars are located in the same part of the HRD. Our 
calculations can be used as a tool to understand the mass 
loss changes of an LBV in terms of changes in Toff during 
such a typical variation (see also Leitherer et al. 1989). 

Observations of LBVs show that for some LBVs that 
undergo typical variations M is increasing from visual 
minimum to maximum, while for others it is the other way 
around: M is decreasing. This "unpredictable" behaviour 
of M during an LBV variation is not a complete surprise, 
if one considers our M values as a function of Toff. We 
have found that in the ranges T^s = 40 000-30 000 K and 
Teff = 20 000-12 500 K, M decreases for a decreasing T^s, 
whereas in the interval between T^s = 30 000 — 20 000 K, 
M increases for a decreasing T^ff . This shows that whether 
one expects an increasing or decreasing M during an LBV 
variation depends on the specific range in Teg between vi- 
sual minimum and maximum. This was already suggested 
by Lamers (1997), albeit a constant value of rj was antic- 
ipated. 

Our present calculations cannot be used to model the 
observed LBV variations, because we have assumed so- 
lar metallicities, whereas the LBVs are known to have an 
enhanced He and N abundance (e.g. Smith et al. 1994). 
Moreover, since most LBVs have already suffered severe 
mass loss in the past, their L*/M* ratio will be higher 
than for normal OB supergiants. This means that LBVs 
are closer to their Eddington limit, which one may expect 
to have an effect on M also. These combined effects ex- 
plain the lack of a consistent behaviour of M for LBV 
variations so far. Especially since it is not sure that 
really remains constant during the variations (see Lamers 
1995). 

8. Summary, Discussion, Conclusions & Future 
work 

We have investigated the nature of the observed jump in 
the ratio Vqo/vcsc of the winds of supergiants near spectral 
type Bl. 

Calculations for wind models of OB supergiants show 
that around T^s = 25 000 K the mass-loss rate M jumps 
due to an increase in the line acceleration of Fe ill below 
the sonic point. This jump in M is found in three different 
series of models. In all cases, the wind efficiency number 
T] = Mvoo/ (L^/c) increases significantly, by about a factor 
of 2 to 3, if Toff decreases from about 27 500 K to about 
22 500 K. Observations show that the ratio Vqo/vosc drops 
by a factor of two around spectral type Bl. Applying these 
values for Voo/vosc, we predict a bi-stability jump in M of 
about a factor of five. So M is expected to increase by 
about this factor between 27 500 and 22 500 K. 

We have argued that the mass loss is determined by 
the radiative acceleration in the subsonic part of the wind, 
i.e. below r ~ 1.03i?*. We found that this radiative accel- 
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eration is dominated by the contribution of the Fe ill Hnes. 
Therefore M is very sensitive to both the metal abundance 
and to the ionization equihbrium of Fe. Our models show 
that the ionization fraction of Fe ill and the subsonic ra- 
diative acceleration increases steeply between Tes= 27 500 
and 25 000 K. This explains the calculated increase in AI 
in this narrow temperature range. 

The exact temperature of the bi-stability jump is 
somewhat ambiguous. Observations indicate that the 
jump occurs around spectral type Bl, corresponding to 
Teff ~ 21 000 K (Lamers et al. 1995). If one would not 
completely trust the value of Vqo/vcsc for the star HD 
109867 (number 91 in Lamers et al. (1995)), because of its 
relatively large error bar, then Tea of the observed jump 
can easily occur at a few kK higher. In fact we cannot 
expect the bi-stability jump to occur at one and the same 
temperature for all luminosity classes, because the jump is 
sensitive to the ionization balance (mainly of Fe III) in the 
subsonic region of the wind and hence to the gravity of the 
star. Our models predict that the jump will occur near T^s 
~ 25 000 K. However, this is sensitive to the assumptions 
of the models: the adopted masses and luminosities and to 
the assumption of the modified nebular approximation for 
the calculation of the ionization equilibrium of iron (see 
Sect.!). 

A more consistent treatment of the ionization and ex- 
citation equilibrium of the Fe-group elements may have 
two effects: i) M predicted from AL may alter, and ii) 
Teff at which the ionization ratio of e.g. Fe iii/iv flips, 
may shift. Nevertheless, in view of the very encouraging 
results using the modified nebular approximation in the 
modeling of UV metal-hne forests (de Koter et al. 1998), 
we expect the error in Tgff at which the dominant ioniza- 
tion of Fe switches from IV to ill to be at most a few kK. 
Furthermore, if a more consistent treatment would yield a 
change in M this would most likely produce a systematic 
shift. Since we are essentially interested in relative shifts 
in M , we do not expect that our conclusions regarding the 
nature of the bi-stability jump would be affected. 

It is relevant to mention that Leitherer et al. (1989) 
calculated atmospheric models for LBVs and suggested 
that the recombination of iron group elements from dou- 
bly to singly ionized stages, which according to them, oc- 
curs around T^s = 10 000 K, can explain M increases 
when LBVs approach their maximum states. We have 
found a Fe iii/iv ionization/recombination effect around 
Teff = 25 000 K for normal supergiants. We also anticipate 
that somewhere, at a lower value of Teff a similar ioniza- 
tion/recombination effect will occur for Fe ii/lii, causing 
a second bi-stability jump. Lamers et al. (1995) already 
mentioned the possible existence of a second bi-stability 
jump around Toff = 10 000 K from their determinations 
of Voo/vcsc, but the observational evidence for this sec- 
ond jump was meagre. Possibly, this second jump is real 
and we anticipate that this second jump could very well 
originate from a Fe ii/iii ionization/recombination effect. 



Furthermore, we have shown that the elements C, N 
and O are important line drivers in the supersonic part of 
the wind, whereas the subsonic part of the wind is domi- 
nated by the line acceleration due to Fe. [| Therefore, we 
do not expect CNO-processing to have a large impact on 
M, but it might have impact on the terminal velocities. 

Finally, we would like to add that our calculations for 
M around Toff — 25 000 K have only been performed 
for one value of L» and H/He abundance. M is ex- 
pected to depend on these stellar parameters, so calculat- 
ing mass-loss rates for a wider range of stellar parameters 
will provide valuable information on the size of the bi- 
stability jump in Vao and M and will allow us to constrain 
its amplitude and location in the HRD. 
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